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Abstract 

Introduction: Due to its remarkable species diversity and micro-endemism, Madagascar has recently been suggested 
to serve as a biogeographic model region. However, hypothesis-based tests of various diversification mechanisms that 
have been proposed for the evolution of the island's micro-endemic lineages are still limited. Here, we test the fit of 
several diversification hypotheses with new data on the broadly distributed genus Eulemur using coalescent-based 
phylogeographic analyses. 

Results: Time-calibrated species tree analyses and population genetic clustering resolved the previously 
polytomic species relationships among eulemurs. The most recent common ancestor of eulemurs was estimated 
to have lived about 4.45 million years ago (mya). Divergence date estimates furthermore suggested a very 
recent diversification among the members of the "brown lemur complex", i.e. former subspecies of £ fulvus, 
during the Pleistocene (0.33-1.43 mya). Phylogeographic model comparisons of past migration rates showed 
significant levels of gene flow between lineages of neighboring river catchments as well as between eastern 
and western populations of the redfronted lemur (£ rufifrons). 

Conclusions: Together, our results are concordant with the centers of endemism hypothesis (Wilme et al. 2006, 
Science 312:1063-1065), highlight the importance of river catchments for the evolution of Madagascar's micro-endemic 
biota, and they underline the usefulness of testing diversification mechanisms using coalescent-based phylogeographic 
methods. 
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Introduction 

Although biodiversity is higher in the tropics than in 
temperate regions, most of our knowledge of species 
dynamics in space and time come from the northern 
hemisphere [1,2]. Climatic changes during the ice ages, 
however, also had profound effects on the history and 
formation of tropical species because cooler and drier 
periods during the Quaternary caused reduction of 
tropical forests and expansion of savannahs [3-5]. As 
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tropical regions are the placeholders and producers of 
great parts of biodiversity, there is an urgent need to 
study those regions [6], and hypothesis-based statistical 
phylogeographic methods are particularly appropriate 
methods for this purpose [7-9]. 

The fourth-largest island of the world, Madagascar, is 
renowned for its exceptional biodiversity and levels of 
endemism [10,11]. New species are still being regularly 
discovered, including plants, reptiles, fishes and mam- 
mals [12]. One hundred percent of amphibians, 90% of 
plants, 92% of reptiles and the primate suborder Lemuri- 
formes are endemic to the island [13], highlighting 
Madagascar's importance for biodiversity studies and con- 
servation efforts [11,14], In addition, a large proportion of 
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Madagascar's extant fauna is micro-endemic to small 
ranges within the landmass of the island [13,15,16]. 

The current understanding of the origin of Madagas- 
car's exceptional faunal biodiversity and endemism is 
that most of the endemic lineages at higher taxonomic 
levels (families and genera) resulted from oversea disper- 
sal from the African or Indian mainland starting about 
65 mya [17], whereas other faunal elements are rem- 
nants of the Gondwanian fragmentation during the Cret- 
aceous when India-Madagascar broke off from Africa 
around 158-160 mya, from Antarctica around 130 mya, 
and the separation of Madagascar from India around 
84-96 mya [16,18,19]. Whereas the origin of these en- 
demic genera and families in Madagascar is well explained 
by irregular colonization events from the African and In- 
dian mainlands, the origin of Madagascar's micro-endemic 



biota is still in debate [15,20]. Several mechanisms have 
been proposed to explain the diversification of Madagas- 
car's extant fauna, recently reviewed by Vences et al. [16]. 

An early model to explain species distributions in 
Madagascar was based on phytogeography, bioclimatic 
zonation of the island and the distribution of lemur spe- 
cies communities [14,21] (Figure lc- d). Following this 
model, the island was separated into eight zoogeographic 
regions and specifically highlighted the importance of 
the western dry and eastern humid habitats, as well as 
the role of some rivers, to further divide similar climatic 
regions [22]. Additional new evidence and changing phy- 
logenies for several taxonomic groups over the last two 
decades, however, revealed considerable discordance be- 
tween these zoogeographic regions and the biogeo- 
graphic separation of Madagascar into an eastern and 
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Figure 1 Maps of Madagascar showing the distribution of Eulemur species, sampling sites and relevant information for different 
diversification hypotheses, a) Distribution of species of the brown lemur complex, formerly considered subspecies of E. fulvus and sampling 
sites. Circles = sampling sites for individuals used in this study, Triangles = Sampling sites of museum specimens, b) Distribution of E. flavifrons, 
E. macaco, E. rubriventer, E.coronatus and £ mongoz. c) Major climatic zones of Madagascar, d) Major eco-geographic regions based on climatic 
zones, e) Illustration of the three highest mountains of Madagascar and associated rivers that are at the base of the centers of endemism (river 
catchments) hypothesis, f) Map showing the centers of endemism (numbered from 1-12) and the retreat dispersal watersheds in between, 
c), d), e) and f) adapted after Vences et al. [16]. 
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western domain [14,22,23], leading to the formulation of 
new hypotheses. 

Wilme et al [15] proposed one hypothesis to explain 
the evolutionary history and regional speciation of Mada- 
gascar's forest biota based on 35,400 occurrence records 
of vertebrate taxa and the watersheds associated with the 
islands rivers. After this so-called centers of endemism 
hypothesis, quaternary paleoclimatic variation played an 
important role in shaping the distribution and speciation 
of the extant Malagasy fauna. Accordingly, during periods 
of glaciation, cooler and dryer climates resulted in more 
arid conditions, forcing forest fragments and forest 
dependent animals to remain isolated along rivers (re- 
fugia). In watersheds with headwaters at low altitudes 
these forest fragments became isolated by intervening 
arid areas, creating centers of endemism, which allowed 
for allopatric speciation and the evolution of micro- 
endemic taxa (Figure If). In contrast, in watersheds of riv- 
ers with sources at high elevation these contracted forest 
fragments could have remained or be newly established at 
higher altitudes and create so-called retreat- dispersal wa- 
tersheds, which would have allowed dispersal along the 
river catchments to neighboring retreat-dispersal wa- 
tersheds (see also [16] for graphical illustration). As 
Madagascar has three major mountains along the east- 
ern highlands above 2000 m (Figure le), and the lar- 
gest rivers of the west (Betsiboka, Tsiribihina and 
Mangoky) as well as of the east (e.g. Manangoro) have 
their headwaters near the summits of those mountains, 
gene flow from the west to the east and vice versa 
would have been possible. 

Pearson and Raxworthy [20] discussed a climatic gra- 
dient model to explain local speciation patterns based 
on current distributions of lemurs, geckos and chame- 
leons, and compared it to the centers of endemism hy- 
pothesis and a biogeographical null model. They found 
concordant distributions with either the centers of en- 
demism or their current climate hypothesis, and sug- 
gested that multiple processes have played a role in the 
diversification of Madagascar's micro-endemic fauna. 

In 2009, Vences et al. reviewed all currently proposed 
diversification hypotheses for Madagascar and formu- 
lated specific predictions to investigate the role of each 
model for the evolution of Madagascar's micro-endemic 
biota. They included five different speciation mecha- 
nisms that are also relevant in other parts of the world, 
which are shortly explained in the following (see Vences 
et al. [16] for details). 

The 'ecogeographic constraint' model is identical to 
the one formulated by Martin ([21], see above) and as- 
sumes that an ecologically tolerant species occurs in dif- 
ferent eco-geographic regions, whereas younger sister 
lineages to the former are more specialized and re- 
stricted to one of the eco-geographic regions (Figure 1 



c-d). Lineages should correspond to eco-geographic re- 
gions and an east- west pattern should be evident. 

A variant of the eco-geographic constraint model, the 
western rainforest refugia' model, assumes that eastern 
species spread into western Madagascar during more 
humid times and became subsequently isolated in rain- 
forest relict areas, which allowed for vicariant speciation. 
No gene flow from west to east can be predicted for this 
mode of speciation. 

The riverine barrier' model assumes rivers to act as 
barriers and allows for allopatric speciation. No gene 
flow between populations or species on both sides of a 
river can be expected from this model, but species on 
one side of the river should be sister species to the ones 
on the other side of the river. 

The montane refugia' hypothesis is based on the as- 
sumption that isolated populations of a widely distrib- 
uted species on high mountains during dry periods later 
diversified due to vicariant divergence. Sister species in a 
phylogeny would be distributed on neighboring massifs 
according to this scenario. 

Finally, the 'river catchments' hypothesis corresponds 
to the centers of endemism hypothesis as proposed by 
Wilme et al. [15]. For species distributed in retreat dis- 
persal watersheds we can expect that gene flow occurred 
several times during pleistocene climatic variations and 
that speciation therefore should have occurred within 
the last ~5 million years [15,16]. As for the 'riverine bar- 
rier hypothesis' species distributed in neighboring retreat 
dispersal watersheds should be sister species in a 
phylogeny. 

Given the various diversification mechanisms, explicit 
hypothesis testing using either the whole Malagasy ecosys- 
tem [16] and/or specific radiations within the extant 
fauna, is now possible (see [7,20,24-27]). The genus of true 
lemurs (Eulemur, [28]) has already been subject to various 
phylogenetic and biogeographic analyses [14,22,29-31]. 
The genus contains 12 species that are distributed over 
the remaining forest fragments of almost the entire island 
of Madagascar (Figure la-b) [32,33]. Seven species, 
namely E. albifrons, E. cinereiceps, E. collaris, E. fulvus, 
E. rufifrons, E. rufus and E. sanfordi, long had unre- 
solved phylogenetic relationships among each other 
and were traditionally classified as subspecies of the 
common brown lemur (E. fulvus) and collectively re- 
ferred to as the 'brown lemur complex' [34], Using 
multiple lines of evidence, Markolf et al. [35] could 
show that all members of the 'brown lemur complex' 
qualify as true species under the general lineage con- 
cept of species [36], supporting an earlier suggestion 
by Groves [37]. 

The species of the 'brown lemur complex' are distrib- 
uted in allopatric populations in a circle-like pattern 
along the remaining forest fragments of the island 
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(Figure la). The only biogeographic zones not inhabited 
by members of the 'brown lemur complex' are the 
central highlands and the south-western spiny forests 
[33]. Eulemur rufifrons and E. fulvus have disjunct 
populations in eastern as well as western parts of the 
island. The remaining members of the genus (Figure lb), 
E. coronatus, E. mongoz, E. rubriventer, E. macaco and 
E. flavifrons occur in sympatry with one of the members 
of the 'brown lemur complex' and exhibit much greater 
genetic divergence among each other and to the members 
of the brown lemur complex [35]. 

Given the broad geographic distribution of eulemurs, 
it is not surprising that the genus Eulemur had an influence 
on the development of several of the above-mentioned hy- 
potheses, including the role of rivers [22,30], the zonation 
into zoogeographic regions [21] or the centers of endemism 
[15]. As the distributions of some species, e.g. E. coronatus, 
E. fulvus, E. sanfordi and E. albifrons [30], are still poorly 
defined, and contemporary distributions do not necessarily 
correspond to distributions during times of speciation, in- 
corporation of phylogeographic approaches, such as gene 
flow models and divergence estimates of species, will help 
to illuminate diversification mechanisms. Thus, the aims of 
this study were two-fold. First, we aimed to resolve the 
phylogeny of the genus Eulemur using multi-locus 
coalescent-based species tree analyses. Second, we wanted 
to infer the predominant speciation mechanisms that 
shaped the evolutionary history of this genus in space and 
time, using coalescent-based phylogeographic methods. 

To this end, we tested the following predictions (see also 
Table 1). For the eco-geographic constraint hypothesis, we 
predicted that distribution of lineages should coincide 
with major Malagasy eco-geographic zones. Furthermore, 
the youngest sister lineage of a group or species should be 
a generalist and occur in different eco-geographic zones, 
We also predicted an east-west phylogeographic pattern 
corresponding to the humid eastern rain forest and the 
western dry forests. 

According to the western refugia hypothesis, we pre- 
dicted no gene flow from west to east. However, this 
model is only relevant for E. fulvus and E. rufifrons, 
which both have populations in the east and the west, as 
well as for E. rufus, which might be a relict population 
of E. rufifrons expanding from the east to the west. 

The riverine barrier hypothesis predicted that sister 
lineages are neighbors and separated by a river. Gene 
flow between sister lineages should be small or absent, if 
rivers are the primary cause of geographic separation 
and divergence. The riverine barrier hypothesis allowed 
specific predictions for all species except E. rubriventer. 

Finally, the river catchment hypothesis predicted that 
lineages occurring in retreat dispersal watersheds are sis- 
ter lineages to lineages in neighboring retreat-dispersal 
watersheds. If Pleistocene glacial maxima and minima 



have been the driving factor for the retreat of popula- 
tions along watersheds, lineages of the brown lemur 
complex must have diverged very recently ( < 5 mya; 
Vences et al. [16]) and watersheds would have allowed 
for gene flow among sister lineages or populations of 
species that occur in eastern as well as in western parts 
of Madagascar, such as E. fulvus and E. rufifrons [15]. 

Results 

Detailed descriptions of the genetic loci used in this 
study are given in Markolf et al. [35]. 

Divergence dates estimation and phylogeny of mtDNA 

Phylogenetic relationships and divergence dates as esti- 
mated from the Bayesian MCMC approach for the 
complete cytb locus are shown in Figure 2 (see also 
Additional file 1: Figure SI + S2). Details about diver- 
gence dates and node support are summarized in Table 2. 
Phylogenetic relationships among higher clades are well 
supported and in agreement with recently published 
phylogenetic relationships among major lineages of the 
Lemuriformes based on multiple genetic loci [38]. Our 
divergence dates, however, are considerably younger for 
deeper nodes than estimated by Perelman et al. [38], but 
correspond to the estimates based on whole mtDNA ge- 
nomes of Finstermeyer et al. [39] that were also used to 
calibrate three of the deeper nodes in the present ana- 
lysis. The most recent common ancestor (MRCA) of all 
eulemurs is estimated to have lived at about 6.15 mya. 
Monophyly is highly supported for the genus Eulemur as 
well as for brown lemur complex (posterior probability 
(pp) = 1) and sister species relationships of E. macaco-E. 
flavifrons (pp = 1) and E. cinereiceps-E. collar is (pp = 1). 
Eulemur rubriventer is the sister lineage to the brown 
lemur complex, but this node is only poorly supported 
(pp = 0.22). The brown lemur complex began to diver- 
sify at the Pliocene-Pleistocene boundary around 1.22- 
3.26 (mean = 2.18) mya. Whereas E. albifrons, E. fulvus 
and E. sanfordi are polyphyletic, the remaining lineages of 
the brown lemur complex, i.e. E. cinereiceps, E. collaris, 
E. rufifrons and E. rufus, are monophyletic for the cytb 
locus (see also [35]). 

Time-calibrated multi locus species tree 

The time-calibrated species tree for the genus Eulemur 
is depicted in Figure 3. Detailed divergence dates and 
posterior probabilities for all clades are given in Table 2. 
Relationships among deeper nodes of the species tree 
correspond to the phylogenetic relationships estimated 
for the cytb locus. Eulemur coronatus, E. macaco and 
E. flavifrons form a sister clade to the remaining eule- 
murs. Eulemur rubriventer is the sister lineage to the 
species of the brown lemur complex. However, this 
node is also not well supported. The monophyly of the 



Table 1 Species-specific predictions (left) and support (right) for different diversification hypotheses modified after Vences et al. [16] for the genus Eulemur 



Species Eco-geographic constraints 



Support Western 
refugia 



Support Riverine barrier 



Support Centers of endemism (river catchments) 



Support 
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General Youngest sister lineage is generalists Partly 
and occurs in different eco-geographic 
zones; Older sister linages are more 
specialized; East (humid)-west (dry) 
phylo-geographic pattern; lineages 
correspond to eco-geographic regions 

Eulemur Lineage distribution corresponds to No 
albifrons eco-geographic regions; east- west 
division 

Eulemur Lineage distribution corresponds to No 
cinereiceps eco-geographic regions; east- west 
division 

Eulemur Lineage distribution corresponds to No 
collaris eco-geographic regions; east- west 
division 

Eulemur Lineage distribution corresponds No 
fulvus to eco-geographic regions; east- 
west division 



Eulemur Lineage distribution corresponds to 
rufifrons eco-geographic regions; east- west 
division 

Eulemur Lineage distribution corresponds to 
rufus eco-geographic regions; east- west 
division 

Eulemur Lineage distribution corresponds to 
sanfordi eco-geographic regions; east- west 
division 

Eulemur Sister lineage to fulvus group 
rubriventer 

Eulemur Lineage distribution corresponds to 
macaco eco-geographic regions; east- west 
division 

Eulemur Lineage distribution corresponds to 
flavifrons eco-geographic regions; east- west 
division 

Eulemur Lineage distribution corresponds to 
coronatus eco-geographic regions; east- west 
division 



No gene flow 
from west to 
east 



No 



Sister lineages occur on either 
side of a river, low gene flow 
between sister lineages 



Sister lineage to E sanfordi 
(Bemarivo), but low gene flow 



Sister lineage to E collaris 
(Mananara) 



Sister lineage to E cinereiceps 
(Mananara) 



Partly 



No 



Yes 



Yes 



Yes 



No 



No 



Yes 



No gene flow 
from west to 
east 

No gene flow 
from west to 
east 

No gene flow 
from E rufus 
to E rufifrons 



Probably no Sister lineage to E rufifrons 
(few data) (Manangoro), E rufus (Betsiboka), 

E albifrons (Mananara), but low 

gene flow 



No 



No 



Sister lineage to E rufus (Tsiribihina), 
but low gene flow 

Sister lineage to E rufifrons 
(Tsiribihina) or fulvus (Betsiboka), 
but low gene flow 

Sister lineage to E albifrons 
(Bemarivo) or fulvus (Mahavavy 
du Nord), but low gene flow 



Sister lineage to E flavifrons 
(Maeverano) 



Sister lineage to E macaco 
(Maeverano) 



Sister species to E macaco + 
E flavifrons (Mahavavy du nord) 



Partly 
Yes 
Yes 
Partly 

Partly 
Partly 
Partly 

Yes 
Yes 
yes 



Lineages occurring in retreat dispersal watersheds Yes 
(RDW) are sister lineages to lineages occurring in 
neighboring RDW; glacial cycles during the 
Pleistocene (< 2.8 mya) allowed gene flow among 
RDW; sister lineages occur in neighboring river 
catchments when their headwaters are at low 
elevations 

Sister lineage to sanfordi {fulvus) + recent Yes 
divergence + gene flow to albifrons and 
fulvus via retreat dispersal watersheds 

Sister lineage to collaris + recent divergence Yes 



Sister lineage to cinereiceps + recent divergence Yes 



Gene flow between east and west 



Gene flow between east and west 



Probably 
(not tested) 



Yes 



Sister linage to rufifrons or fulvus + recent Yes 
divergence + gene flow to fulvus and/or rufus 

Sister lineage to albifrons (fulvus) + recent Yes 
divergence + gene flow to albifrons via 
retreat dispersal watersheds 



Sister lineage to E flavifrons + recent divergence Yes 



Sister lineage to E macaco + recent divergence Yes 



Sister lineage to macaco/flavifrons + recent Yes 
divergence 



The western refugia hypothesis allowed only to formulate predictions for some of the species. RDW: Retreat dispersal watershed (after Wilme et al. 



[15]). 
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Millions of years before present 

Figure 2 Simplified Bayesian tree of 121 Eulemur individuals and seven outgroups of the cytochrome B with divergence date estimates 
for well supported nodes. The mean age is given in millions of years and 95% credibility intervals are indicated by blue bars. A geological time 
scale is shown at the top. Details of divergence date estimates are given in Table 2. 



Table 2 Bayesian divergence date estimates in million of years 



Node 




Cytochrome B 






Species tree 




Mean 


95% HPD 


Prob 


Mean 


95% HPD 


Prob 


Chiromyiformes + Lemuriformes- Lorisiformes* 


56.71 


51.2- 62.34 


1 








Chiromyifromes - Lemuriformes* 


46.77 


39.77- 53.84 


1 








Propithecus - Lemuridae + Cheirogaleidae* 


27.68 


22.54- 33.21 


1 








Lemuridae - Cheirogaleidae 


22.34 


14.88- 28.95 


0.85 








Lemuridae 


14.56 


10.92- 22.76 


0.84 








Lemur cotta- Hopolemur griseus 


9.31 


14.35- 14.27 


0.89 








MRCA Eulemur 


6.15 


3.6- 8.89 


1 


4.45 


3.26- 5.68 


1 


MRCA £ coronotus + E. macaco + E. flavifrons 


4.46 


2.42- 6.8 


0.87 


3.84 


2.65- 5.05 


0.58 


MRCA E. macaco + E. flavifrons 


2.04 


0.91-3.31 


1 


1.15 


0.6- 1.71 


1 


MRCA fulvus group + E. rubriventer + E. mongoz 


4.55 


2.61-6.81 


0.96 


2.86 


1.83- 3.91 


1 


MRCA fulvus group + E. rubriventer 


4.06 




0.22 


2.24 


1.16- 3.32 


0.6 


MRCA fulvus group 


2.18 


1 .22- 3.26 


1 


0.93 


0.33- 1 .43 


0.98 


MRCA E. albifrons, E. fulvus, E. sanfordi, E. rufifrons, E. rufus 








0.35 


0.22- 0.51 


0.9 


MRCA E. cinereiceps + E. collaris 


0.8 


0.3-1.38 


1 


0.51 


0.22- 0.79 


0.91 


MRCA E. rufifrons + E. rufus 








0.17 


0.08- 0.28 


0.98 


MRCA E. fulvus + E. albifrons + E. sanfordi 








0.27 


0.19- 0.36 


0.86 


MRCA E. albifrons + E. sanfordi 








0.09 


0.04- 0.14 


1 



The mean, 95% credibility intervals (95% HDP) and node supports (Prob) are given for the analyses of the cytochrome B and the species tree estimation from 
multiple loci. Missing values (-) are due to taxa that were not included in the species tree estimation, low support or discordance among the gene tree of the 
cytochrome B and nodes estimated from the combined analysis of multiple loci. MRCA = Most Recent Common Ancestor, * = time calibrated nodes from Table 6. 
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Millions of years before present 

Figure 3 Time-calibrated species tree of the genus Eulemur based on two mitochondrial and three nuclear loci. Black solid lines show a 
single combined tree estimated from 80.004 species trees. Numbers depict posterior probabilities of each node. Gene tree species tree 
discordance is illustrated by 10.000 colored trees of the posterior distribution in the background. Blue: Most popular topologies, Yellow: 2nd most 
popular topologies, Green: 3rd most popular topologies. A geological time scale is given at the top. Details of species divergence dates are given 
in Table 2. 



brown lemur complex is well supported (pp = 1) as are 
the sister group relationships of E. collaris and E. cinereiceps 
(pp = 0.91), E. rufus and E. rufifrons (pp = 0.96), and 
E. albifrons and E. sanfordi (pp = 1). The sister group 
relationship of E. fulvus to E. albifrons and E. sanfordi 
is supported by a posterior probability of pp = 0.86. 
Species divergence dates are similar but slightly younger 
compared to the cytb locus and 95% credibility intervals 
are smaller for the multi locus analysis. The most recent 
common ancestor of the genus Eulemur was estimated to 
have lived at about 4.45 (3.26-5.68) mya. Eulemur macaco 
and E. flavifrons diverged from E. coronatus about 3.84 
(2.65-5.05) mya. Eulemur macaco and E. flavifrons di- 
verged about 1.15 (0.6-1.71) mya. Eulemur mongoz di- 
verged from E. rubriventer and the members of the 
brown lemur complex about 2.86 (1.83-3.91) mya. The 
split between E. rubriventer and the members of the 
brown lemur complex was dated at 2.24 (1.16- 3.32) 
mya. The MRCA of the brown lemur complex was es- 
timated at 0.93 (0.33-1.43) mya. The clade was then 
split into the two most southern species, E. cinereiceps 
and E. collaris that diverged 0.51 (0.22-0.79) mya, and 
the remaining species of the brown lemur complex 



that diverged 0.35 (0.22 0.51) mya into two groups, 
one containing E. rufus and E. rufifrons and one con- 
taining E. albifrons, E. fulvus and E. sanfordi. Splits of 
E. rufus-E. rufifrons and E. albifrons-E. sanfordi were 
estimated at only 0.17 (0.08-0.28) mya and 0.09 (0.04 
0.14) mya, respectively. Diversification of the brown 
lemur complex occurred during the last -1.5 million 
years of the late Pleistocene. The species tree estimated 
without the PAST fragment resulted in similar diver- 
gence date estimates and similar phylogenetic relation- 
ships among most of the clades (see Additional file 1: 
Figure S2). However, the positions of E cinereiceps, E. collaris 
and E. fulvus were different, and posterior probabilities for 
all clades are considerably lower. 

Nuclear genetic population structure 

Genetic population structure of three nuclear loci of the 
members of the brown lemur complex as estimated with 
STRUCTURE and DAPC in Markolf et al. [35] plotted 
on a map of Madagascar is depicted in Figure 4. For the 
STRUCTURE results of K = 3 populations, individuals 
from the east cluster with individuals from the west, and 
a clear south to north structure is evident. Assignment 
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Figure 4 Geographic plot of nuclear genetic population structure of species of the brown lemur complex as inferred by Markolf et al. [24] 
using STRUCTURE (K = 3) and Discriminant Analysis of Principal Components (DAPC). Pies represent individuals. Colors represent assignment 
probabilities of individuals to populations (STRUCTURE, left) or species (DAPC, right). Species colors for the DAPC analysis are given in the color legend. 
Please note that the color legend is only relevant for the map on the right. Pies correspond only roughly to the sampling locality. 



probabilities of the DAPC supports the sister group rela- 
tionship of E. sanfordi and E. albifrons as estimated in 
our species tree in northern Madagascar as well as sig- 
nificant differentiation of nuclear genes of E. fulvus and 
E. rufus. Western and eastern populations of E. rufifrons 
show mixed nuclear genetic composition. Eulemur col- 
laris individuals in the southeast are best separated from 
the others based on nuclear genetic data although some 
admixture exists with eastern E. rufifrons. 

Model-based phylogeography 

Marginal likelihoods corresponding to Bayes factors and 
relative model probabilities of the different migration 
models for three population/species combinations are 
reported in Table 3,4,5. In all cases coalescent simula- 
tions favored the more complex model of a full migra- 
tion matrix between populations/species over more 
simpler models of panmixia, uni-directional or no gene 
flow models. Although we tested all possible combina- 
tions for the dyads or triads, we only report the models 
that had biological relevance in terms of the potential 
speciation mechanisms mentioned above. Past immigra- 
tion rates were high, especially for the migration model 
of eastern and western populations of E. rufifrons. How- 
ever, as we did not aim to interpret and assess the exact 
number of migrants or the effective population sizes, 
demographic parameters of 0 and M over all loci for the 
best models are reported in Additional file 2: Table S2. 



Here, our aim was to test the prediction of past gene flow 
between sister lineages of the species tree or species that 
occur in disjunct populations on both sides of the island. 
All three models clearly rejected panmixia or the no gene 
flow models (p < 0.001) and favored a full migration 
matrix model (gene flow in all directions) with a relative 
probability to all other models of 0.999. 

Results for the specific predictions for different diversi- 
fication hypotheses are summarized in Table 1. The 
combination of species divergence dates, which corres- 
pond well to the climatic variations during glacial cycles 
in the late Pleistocene, sister group relationships as esti- 
mated from the species tree, and Bayes Factor compari- 
sons of gene flow models are highly concordant with the 



Table 3 Log marginal likelihoods (ImL) and log Bayes 
factor (LBF) comparisons for different migration models 
for western and eastern populations of E. rufifrons 



Model 


BA ImL 


LBF 


Model prob 


Model rank 


full migration matrix 


-3056.85 


0 


1 


1 


panmixia 


-3129.01 


-72.16 


<0.001 


4 


no gene flow 


-3193.61 


-136.76 


<0.001 


5 


west to east 


-3085.35 


-28.5 


<0.001 


3 


east to west 


-3084.74 


-27.89 


<0.001 


2 



The log marginal likelihood is given as Bezier approximation score (BA ImL). 
LBF shows differences between the best and all other models. The model 
probability (Model prob) shows the probability of each model being the 
correct model relative to the others. 
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Table 4 Log marginal likelihoods (ImL) and log Bayes 
factor (LBF) comparisons for different migration models 
for E. fulvus, E. rufifrons and E. fulvus 



Model 


BA ImL 


LBF 


Model prob 


Model rank 


full migration matrix 


-4786.19 


0 


1 


1 


panmixia 


-5032.93 


-246.74 


<0.001 


2 


no gene flow 


-5137.37 


-351.18 


<0.001 


4 


rufifronsofulvus 


-5190.23 


-404.04 


<0.001 


5 


fulvusorufus 


-5084.62 


-298.43 


<0.001 


3 


rufifronsorufus 


-5227.18 


-440.99 


<0.001 


6 



Table 5 Log marginal likelihoods (ImL) and log Bayes 
factor (LBF) comparisons for different migration models 
for E. albifrons, E. fulvus and E. sanfordi 



Model BA ImL LBF Model Model 

rank 

prob 



full migration matrix 


-4278.23 


0 


1 


1 


panmixia 


-4498.92 


-220.69 


<0.001 


3 


panmixia albifrons/ 
sanfordi 


-4403.69 


-125.46 


<0.001 


2 


no gene flow 


-4887.55 


-609.32 


<0.001 


5 


E. albifrons <> E. 
sanfordi 


-4518.64 


-240.41 


<0.001 


4 



center of endemism hypothesis. In contrast, we found no 
or only limited support for any of the other hypotheses. 

Discussion 

In this study we explored the evolutionary history of the 
genus Eulemur in space and time and could resolve the 
previously polytomic phylogenetic relationships among 
members of the group. Divergence date estimates indi- 
cate that the MRCA of the genus Eulemur is estimated 
to have lived -4.45 mya and that diversification among 
the members of the fulvus group happened during the 
Pleistocene. Additional comparisons of gene flow models 
among sister lineages favored full migration models over 
panmixia, uni-directional or no gene flow models. After 
discussing the validity of our phylogeographic analyses 
we will discuss the fit of our data to the different diversi- 
fication hypotheses proposed for the evolution of micro- 
endemsim in Madagascar as well as the suitability of our 
approach to other radiations endemic to the island. 

Phylogeography of eulemurs 

The present analyses clearly suggest a Pleistocene origin 
for members of the brown lemur complex as well as for 
E. macaco and E. flavifrons. Divergence dates estimated 
for the cytb locus were slightly older than divergence 
dates for the species tree analysis. This can be explained 



by the smaller effective population size of mtDNA 
compared to nuclear DNA [40] and the fact that gene 
divergence will occur prior to species divergence, and 
divergence dates estimated from single gene trees will 
necessarily overestimate divergence times [41,42]. As 
time-calibrated species trees provide more realistic es- 
timates of species divergence [43] our divergence date 
estimates provide a more realistic picture than previous 
analyses based on single genes or concatenated genes. 

As there are no fossil calibrations points available for 
lemurs [32,44], we used calibrations points from a recent 
study based on complete mitochondrial genomes [39] to 
calibrate our tree for the cytb locus and used the esti- 
mated clock rate from this analysis for the calibration of 
the species tree. As calibration points in [39] were based 
on several dated primate fossils, the clock rate was 
allowed to vary among the remaining loci and the ap- 
plied substitution rate of 0.0138 substitutions/per site/ 
per million years is close to the 2% evolutionary rate for 
vertebrate mtDNA [45], the present divergence date esti- 
mates should therefore not be dramatically over- or 
underestimated. Although accuracy of molecular diver- 
gence dates should not be taken as obsolete, because di- 
vergence date estimations are particularly difficult for 
lemurs due to branch rate variation and the lack of 
lemur fossils [46], a very recent divergence of the brown 
lemur complex in the last four million years is in agree- 
ment with other recently published studies, (e.g. 3.47 
mya (2.58- 4.40) in [39]; 3.1 mya (2.77- 4.04) in [44]; 
2.91 mya (1.57- 4.27) in [38]. 

Simulation studies revealed that three loci combined 
with multiple gene copies per lineage are sufficient to re- 
solve a species tree with high accuracy even of recently 
diverged radiations [43,47-49]. Furthermore, it has been 
shown that concatenation of different genes can lead to 
substantial errors in phylogeny estimation [50]. Although 
the number of gene copies per lineage varied considerably 
between lineages (see Additional file 3: Table SI) because 
our sampling was focused on the members of the brown 
lemur complex, the present species tree analysis using five 
loci represents the most complete phylogeny for the genus 
Eulemur so far. Posterior probability values ranged from 
0.86-1.00 for the phylogenetic relationships among the 
species of the brown lemur complex, which could not be 
resolved in previous studies based on single gene trees or 
concatenated genes [22,32,44,51]. Inclusion of the PAST 
fragment without linking the tree partitions, as suggested 
for mitochondrial DNA in species tree analyses, did not 
introduce any bias to the present phylogeny. As depicted 
in Additional file 1: Figure SI, the phylogenetic relation- 
ships of the PAST fragment are completely concordant 
with the phylogenetic relationships estimated for the cytb 
locus. Eulemur albifrons, E. fulvus and E. sanfordi are 
polyphyletic and E. rufus is a sister group to a clade 
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consisting of E. rufifrons, E. albifrons, E. fulvus and E. san- 
fordi Exclusion of the PAST fragment in species tree ana- 
lysis, however, resulted in a different topology, but 
consistent pattern for deeper nodes. Although both mito- 
chondrial genes did neither find a sister group relationship 
between E. albifrons and E. sanfordi nor between E. rufi- 
frons and E. rufus, the inclusion of three nuclear loci 
seems to support the close relationships among these 
taxa. A sister group relationship of E. albifrons and E. 
sanfordi is also supported by Bayesian nuclear struc- 
ture analysis for K = 3 as shown in Figure 4. 

Bayes factor comparisons of coalescent simulations for 
different phylogeographic models among sister groups 
left little room for misinterpretations of the prevailing 
migration pattern. All three model comparisons consist- 
ently rejected panmixia and the no gene flow model in 
favor of a full migration model among lineages. This is 
highly consistent with several events of gene flow be- 
tween members of adjacent retreat dispersal watersheds 
and the centers of endemism hypothesis [15]. Rejection 
of panmixia furthermore supports the delimitation of 
the members of the brown lemur complex as distinct 
species, as suggested recently by Markolf et al. [24], 
despite a high degree of incomplete lineage sorting due 
to past migration events among lineages during the 
Pleistocene. 

Eco-geographic factors 

The eco-geographic constraints hypothesis can be rejected 
as a general model for the diversification of the genus 
Eulemur. Only three species, E. coronatus, E. rufus and 
E. sanfordi are exclusively distributed in one of the 
eco-geographic zones (Figure 1). However, the position 
of E. rubriventer as the sister lineage to all species of 
the brown lemur complex, and the fact that E. rubriventer 
is distributed along the entire east coast, suggest the possi- 
bility that ecological factors also played a role during the 
initial diversification of the brown lemur complex. If the 
phylogenetic position of E. rubriventer is correct, one 
could hypothesize that populations of the much more 
broadly distributed E. rubriventer had to retreat to isolated 
mountain refugia during cooler and drier periods. Individ- 
uals adapted to more arid conditions, however, could have 
descended from mountain refugia to lower elevations, 
forming the MRCA of the members of brown lemur com- 
plex. This is highly speculative, but is supported by the 
fact that E. rubriventer is normally found at higher eleva- 
tions than sympatric species of the brown lemur complex 
(Markolf, pers. observation, [34]). However, with the 
data currently at hand, this notion is impossible to test, 
not the least because the position of E. rubriventer was 
also one of the least supported in the present phyl- 
ogeny. Although we did not include any ecological var- 
iables in the present analysis, the adaption of E. fulvus 



and E. rufifrons to eastern and western regions with 
very different climatic conditions does not support the 
model of ecogeographic constraints as a general model 
for Eulemur diversification. 

Western refugia 

The western refugia hypothesis predicted no gene flow 
from western to eastern populations. In the present 
dataset, this hypothesis was only biologically relevant for 
western and eastern populations of E. rufifrons, E. fulvus 
and E. rufus, which could potentially be a western relict 
population of eastern E. rufifrons. However, the gene 
flow models clearly reject the predictions of no gene 
flow from west to east for E. rufifrons and E. rufus. Un- 
fortunately, we could not test gene flow between eastern 
and western populations of E. fulvus, because we had 
only two geographically disjunct individuals from the 
west. However the nuclear genetic structure results and 
the phylogeny of the cytb locus (see also [35]) suggested 
gene flow between east and west also for E. fulvus. 

Riverine barriers 

The riverine barrier hypothesis predicted sister lineages 
on either side of a river. This pattern is true for all eule- 
murs based on our genetic sampling and the species tree 
except for E. rubriventer. However, the amount of gene 
flow between sister species that occur on both sides of 
the river is not concordant with a hypothesis that pre- 
dicts rivers as the primary force for the physical separ- 
ation of species. Furthermore, there is evidence that 
large rivers do not form a barrier for several species. 
Eulemur mongoz, for example, is distributed on both 
sides of the Betsiboka, the largest river of Madagascar. 
Goodman and Ganzhorn [30] evaluated the role of rivers 
and the distribution of eulemurs in the eastern rainforest 
and also found no support for the riverine barrier hy- 
pothesis based on eulemur distributions for most taxa. 
Eulemur albifrons and E. fulvus, for example, do not 
have a riverine barrier and might occur in parapatric or 
sympatric populations [30,34], and E. fulvus occurs 
south of its supposed riverine barrier, the Manangoro 
[52]. Therefore, it is highly unlikely that the riverine bar- 
rier hypothesis can explain the diversification and present 
distribution of the genus Eulemur alone. 

Centers of endemism 

Our data broadly support the centers of endemism hy- 
pothesis as the main force in driving Eulemur diversity. 
The prediction of sister species relationships among 
neighboring retreat-dispersal watersheds could be con- 
firmed with high support for all higher nodes in the 
Eulemur phylogeny. Furthermore, the timing of speci- 
ation is concordant with the time of climatic variations 
during glacial cycles of the Pleistocene. As retreat and 
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dispersal to refugia at higher elevations would have hap- 
pened several times during the Pleistocene [15], high 
levels of gene flow among sister species occurring in 
neighboring retreat-dispersal watersheds can be expected 
and were confirmed by our phylogeographic models. Eule- 
mur rubriventer is again the only taxon that shows no 
concordance whatsoever with river catchment hypothesis. 
Fine scale genetic sampling of E. rubriventer along its dis- 
tribution would be necessary to test whether mountain re- 
fugia shaped the demographic history of this species. 

The lack of concordance of E. rubriventer with the 
center of endemism hypothesis also highlights an unrealis- 
tic assumption that one speciation mechanism or diversifi- 
cation hypothesis can and must explain the diversification 
pattern of an entire genus or all radiations endemic to 
Madagascar. Although it might be less important for the 
diversification of the genus Eulemur, the montane refugia 
hypothesis, for example, could be shown to explain pat- 
terns of species richness and endemism in cophyline frogs 
[27]. Furthermore, climatic gradients had probably im- 
portant influences on the diversification of several chame- 
leons, geckos and also lemurs [20]. 

Testing diversification mechanisms with unknown ancestral 
distributions 

It has been shown repeatedly in all major primate radia- 
tions that climatic fluctuations during the Quaternary 
had a fundamental influence on the diversification of 
several primate genera [53-58]. This study, however, rep- 
resents the first example of explicit hypothesis-based 
testing of the diversification mechanism of an endemic 
primate radiation. Our approach using coalescent simu- 
lations was particularly useful because exact distribu- 
tions of Eulemur species are still poorly defined and 
todays distribution must not necessarily correspond to 
the distribution of lineages during speciation events. 
Our geographically broad-scale genetic sampling, how- 
ever, should compensate for uncertainty of ancestral 
lineage distributions. Eulemur sanfordi's distribution, for 
example, is supposed to be restricted to the centers of 
endemism 1 and 12 of Wilme et al. [15] (Ankarana and 
Vohimarina after (Wilme et al.[59])) with the Mana- 
mbato river as its southern limit [34]. However, it can be 
assumed that E. sanfordi had a much wider distribution 
in the past. Evidence comes from a museum sample col- 
lected south of the Manambato close to Vohemar that 
corresponds phenotypically to E. sanfordi and clusters 
with E. sanfordi/E. albifrons based on mitochondrial 
DNA [35] as well as a sample (ID = 491, Additional file 3: 
Table SI) north of the Bemarivo, which is more likely to 
be E. sanfordi based on nuclear genetic assignment prob- 
ability [35]. Unfortunately, we do not have phenotypic in- 
formation on this individual. Additionally, E. coronatus, 
which occurs in sympatry with E. sanfordi, and is also 



supposed to have its southern distributional limit at 
the Manambato river [34], was found at the same locality 
(Anjombalava, samples 490 and 494, Additional file 3: 
Table SI) north of the Bemarivo. We can therefore assume 
that the distribution of E. sanfordi was extended to adja- 
cent RDWs Mahavavy and Bemarivo [59], which allowed 
gene flow to neighboring RDWs during the Pleistocene. 
Our data clearly favored a gene flow model over a model 
of panmixia of E. albifrons and E. sanfordi or a model of 
complete isolation of the latter two, illustrating the power 
of molecular coalescent-based approaches despite un- 
known ancestral distribution to test phylogeographic 
hypotheses. 

Methods to test phylogeographic hypotheses are diver- 
sifying rapidly [7,8], and we are aware of the fact that 
there are several methods, e.g. ecological niche modeling 
approaches [60,61], approximate Bayesian computations 
(ABC) [62,63] or isolation with migration models (IMa) 
[64], that could be additionally applied to the present 
data set to further explore the evolutionary history of 
this group. However, time-calibrated species tree ana- 
lyses and Bayes factor comparisons of gene flow models 
as applied here, using several different model compari- 
sons, could clearly answer our questions concerning di- 
versification of the genus Eulemur in space and time and 
had the advantage over other methods in reducing the 
amount of demographic parameters that have to be esti- 
mated in parallel from the data, especially when the 
number of species is high and computational effort 
would be immense [65]. 

Madagascar as a biogeographic model region 

As previously suggested [16], time is overdue to use 
Madagascar as a biogeographic model region, and to 
conduct hypotheses-based testing of phylogeographic 
pattern among the many endemic lineages to infer speci- 
ation mechanisms that shaped this islands stunning bio- 
diversity. Madagascar is particularly suitable as a model 
region of species diversification because data and sam- 
ples can be collected within the borders of one country, 
which has practical advantages concerning the adminis- 
trative procedures necessary to sample and export bio- 
logical material of CITES listed taxa [16]. Furthermore, 
its high species richness and endemism, together with a 
relatively simple geographic structure of the island, but 
pronounced climatic variations from east to west, to- 
gether with pronounced regional ecotones allows to test 
recurring patterns in several different animal and plant 
radiations in a relatively small geographical area isolated 
from other continental landmasses for a long time. As 
different taxa diversify at different times, several diversi- 
fication mechanisms may have influenced even single ra- 
diations as was also evident from our analysis. 
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Our approach, however, could be easily adapted to 
other endemic radiations of the island that have been 
less involved in the initial formulation of different bio- 
geographic models for Madagascar. It would be particu- 
larly interesting for species that have more restricted 
distributions than the Eulemur species. Genetic data 
already exists for various lineages and genomic resources 
for non-model organisms are increasing rapidly [46,66]. 
Sister lineages of mouse lemurs, for example, showed 
considerable correspondence with the centers of endem- 
ism hypothesis [67], however we do not yet know the 
time of species divergences and, whether they corres- 
pond to major climatic events during the Pleistocene. 
Although the accuracy of species trees, for example, de- 
pends on a optimal range of the number of loci, individ- 
uals and sequence length [68], phylogeographic studies 
can also test diversification hypotheses on a smaller geo- 
graphical scale, as recently shown for northern popula- 
tions of Daubentonia madagascariensis [69] or frogs of 
the genus Mantella [70]. The application of hypothesis- 
based tests on speciation mechanisms to more single 
Malagasy radiations in the future will allow to infer the 
"global" patterns of diversification of Madagascar's bio- 
diversity by integrating multi-locus phylogenies, eco- 
logical niche modeling and GIS approaches in a 
comparative framework [7]. This in turn could help to 
understand the many ways that have shaped biological 
diversity in other regions of the planet. The future of 
phylogeography seems promising due to the advances 
in sequencing technology and statistical modeling 
techniques [8]. However, investigating mechanisms of 
species diversification needs case-specific formulations of 
predictions, which can then be tested with coalescent- 
based phylogeographic techniques [9,71] and/or GIS mod- 
eling techniques [7,60]. 

Conclusions 

We conclude that the diversification of the genus Eulemur 
was shaped by climatic variation during the Pleistocene, as 
suggested by the centers of endemism hypothesis [15]. This 
result highlights the importance of river catchments for the 
evolution of Madagascar's large number of microendemic 
lineages. Nevertheless, other diversification mechanisms, 
such as the role of montane refugia, local or regional cli- 
matic variations or a combination of several different 
forces should not be neglected and could well have played 
a role in the diversification of other radiations on the is- 
land. However, testing these models with genetic data re- 
quires a priori formulated predictions as well as a dense 
sampling design for the lineages under investigation. 

Material and methods 

Genetic data of wild populations of eulemurs collected by 
Markolf et al. [35] and Pastorini et al. [22] were used to 



estimate divergence times and phylogenetic relationships 
for single gene trees as well as for a multi-locus species 
tree. Details of DNA extraction and sequencing have been 
described in detail elsewhere [35]. Nuclear population 
structure of the brown lemur complex as estimated in 
Markolf et al. [35] was plotted on a map of Madagascar 
and gene flow models were compared using a Bayesian 
approach as implemented in migrate-n [72] . 

Divergence date estimation and mtDNA phytogeny 

Sequence data of the complete cytochrome b (1140 bp) 
of 121 Eulemur individuals were used to simultaneously 
estimate phylogeny and divergence times in a Bayesian 
MCMC approach using a relaxed molecular clock as im- 
plemented in Beast version 1.7.5 [73]. Seven additional 
outgroup taxa were included in the analysis. As there 
are no fossil calibration points available for lemurs 
[32,44], calibrations were based on molecular evidence 
from a phylogeny of complete mitochondrial genomes of 
primates [39] as depicted in Table 6. A HKY + I + G 
substitution model was chosen as suggested by Akaike s 
Information Criterion of JModeltest v2 [74]. A birth- 
death process and an uncorrelated log-normal relaxed 
clock with a broad normal prior distribution for the 
mean of the branch rates (ulcd.mean = 0 - <*>) was as- 
sumed. Fifty million generations were run with param- 
eter sampling at every 5.000 generation resulting in 
10.001 trees. 

The adequacy of the burn-in was assessed by visual in- 
spection of the trace of the parameters using Tracer 
v.1.5 [75]. Tree Annotator vl.7.5 was used to discard 
2.500 trees as burn-in and to calculate a maximum clade 
credibility tree of the remaining 7.501 trees. 

Time calibrated multi-locus species tree 

The multi-species coalescent approach implemented in 
*BEAST vl.7.5 was used to infer a species tree for the 
genus Eulemur based on one mitochondrial, three nu- 
clear loci published by Markolf et al. [35] and one mito- 
chondrial locus published by Pastorini et al. [22]. The 
numbers of sequences included were 109 for the cytb 
locus, 147 for the eno locus, 125 for the vwf locus, 120 
for the nramp locus and 53 for the past fragment, result- 
ing in a total number of 554 sequences. Both alleles were 



Table 6 Calibrated nodes, means, standard deviation (SD) 
and 95% confidence intervals in million of years used for 
divergence date estimates of the cytochrome b tree 



Calibration node 


Mean +/- 


- SD 


95% range 


Chiromyiformes + Lemuriformes- 
Lorisiformes 


57.09 +/- 


- 4.2 


50.18- 64 


Chiromyiformes - Lemuriformes 


47.38 +/- 


3.99 


40.82- 53.94 


Propithecus- Lemuridoe 


27.76 +/- 


- 3.1 


22.66- 32.86 
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used for all nuclear loci. * BEAST simultaneously estimates 
gene trees and species trees under the multi species co- 
alescent [47]. As the model assumes that discordance of 
gene trees is based solely on incomplete lineage sorting, 
we had to exclude potential and known hybrids prior to 
analysis (see Additional file 3: Table SI). Potential hybrids 
were determined via discriminant analysis of principal 
components (DAPC) [35]. Exclusion of individuals re- 
sulted in incomplete taxon sampling for some of the loci 
for E. cinereiceps and E. flavifrons. As *BEAST requires at 
least one sequences per species per locus, we included the 
2.400 bp (PAST) fragment of mtDNA published by Pastor- 
ini et al. [22,76] to have sufficient genetic information for E. 
cinereiceps and E. flavifrons. Dummy sequences ( ? = un- 
known state) were coded for the nramp and vwf loci for E. 
cinereiceps and for all three nuclear loci for E. flavifrons. 
Tree, substitution and clock models were unlinked for all 
partitions. As tree partitions of two mitochondrial genes 
should be linked in *BEAST analyses, because mtDNA 
lacks recombination among genes, we calculated two separ- 
ate species trees, once with and once without the PAST 
fragment. Linking tree partitions for the two mtDNA genes 
was not possible, because sample sizes of the cytochrome B 
of Markolf et al. [35] and Pastorini et al. [22] were too 
different. 

To calibrate the species tree in units of million of 
years, we set the clock rate of the cytb locus to the esti- 
mated substitution rate (0.0138) as revealed by the previ- 
ous divergence time analysis of the cytb locus. The clock 
rates of the other loci were allowed to vary. The analyses 
were run with a Birth-Death prior and substitution models 
as indicated by jModeltest v2 (cytb = HKY + I + G, eno + 
vwf = GTR + I, nramp = HKY + G, PAST = GTR + G). 
For both analyses, we ran four separate runs of 30 million 
generations each and a sampling of parameters every 
1.000 generation, resulting in 30.001 trees for each run. 
Convergence of the MCMC runs, adequacy of the burn-in 
and effective sample size (ESS >200) were assessed using 
the combined log.files in Tracer v. 1.5. Trees of separate 
runs were combined using LogCombiner v. 1.7.5 discard- 
ing one third (10.000) of the trees as burn-in for each run. 
Trees of the four separate runs were combined using 
LogCombiner, and TreeAnnotator was used to calcu- 
late the final species tree from 80.004 trees. DensiTree 
[77] was additionally used to visualize gene tree species 
tree discordance using 10.000 trees from the posterior 
distribution. 

Geographical visualization of nuclear population structure 

Nuclear genetic population structure of the members of 
the brown lemur complex estimated in Markolf et al. 
[35] based on a genotype matrix of three nuclear genetic 
loci was plotted on a map of Madagascar, using the on- 
line platform PhyloGeoViz [78]. PhyloGeoViz was 



originally designed to plot haplotype or allele frequencies 
as proportions of pies on a map. However, geo-referenced 
pie charts can also be constructed using assignment prob- 
abilities of individuals to populations inferred from genetic 
clustering methods such as STRUCTURE [79] or Discrim- 
inant Analysis on Principal Components (DAPC) [80]. In- 
dividual assignment probabilities of STRUCTURE for K = 
3 and DAPC (see [35]) were plotted separately on a map 
of Madagascar to geographically visualize nuclear genetic 
population structure. Due to the difficulties of visualizing 
multiple individuals from the same location, the geo- 
graphic positions of pie charts correspond only roughly 
with the sampling site. 

Model-based phylogeography 

Log marginal likelihood comparisons (Bayes factors) of 
coalescent simulations were used to assess the fit of the 
data to different phylogeographic models following the 
approach of Beerli & Palczewski [65] implemented in 
the software MIGRATE-n v3.5.1 [72]. Three different 
model comparisons were conducted following the spe- 
cies tree relationships among eulemurs. Eulemur collaris 
and E. cinereiceps were not included because of small 
sample size. Model comparisons were conducted be- 
tween western and eastern populations of E. rufifrons, 
between E. fulvus, E. rufifrons and E. rufus and finally be- 
tween the three northern species of E. fulvus, E. albifrons 
and E. sanfordi. The three nuclear genetic loci and the 
complete cytb locus of Markolf et al. [35] were used for 
the analyses. The mutation rates for the three nuclear loci 
were scaled to 0.25, comparable to mtDNA, using the in- 
heritance scalar in MIGRATE-n to allow for easy inter- 
pretation of multi-locus parameters. Markers were run 
with a F84 substitution model and transition/transversion 
ratios of 13.1 (cytB), 2.3 (eno), 2.3 (nramp) and 3.1(vwf) as 
indicated by jModeltest v2. Mutation rate was set to con- 
stant, as suggested for most analyses by the user manual 
of migrate-n [72]. Bayesian analysis consisted of one long 
chain with 10.000 recorded parameter steps, a sampling 
interval of 100 and a burn-in of 250.000 (25 %). We used 
Metropolis Hastings sampling and eight statically heated 
chains at their default temperatures simultaneously in 
each run to effectively explore the parameter space. Uni- 
form prior distributions for 0 and M were assumed. 

To compare models, scaled log Bayes factors were cal- 
culated by subtracting the highest value of the log mar- 
ginal likelihoods (lmL) (Bezier curve approximation) 
from lmL values of each model. The probability of the 
model in relation to all other models tested was then 
calculated by dividing the Bayes factor by the sum of 
all Bayes factor scores from all models following Kass 
& Raftery [81]. For all three model combinations, we 
tested all possible combinations. However, we report 
and describe only those that are biological meaningful 
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in terms of the species distribution and the island 
geography. Those were a full migration matrix model 
(gene flow in all directions among all populations), a 
panmixia model, where populations are treated as one 
panmictic population, and a no gene flow model by setting 
M to a constant value of 0.1 migrant per generation. For 
eastern and western populations of E. rufifrons, we add- 
itionally included a model with asymmetrical gene flow 
between east and west. For the three species comparison 
of E. rufifrons, E. fulvus and E. rufifrons, we additionally 
included models that predict only gene flow between two 
of these populations, which could be equally likely to a full 
migration matrix model based on the distribution of the 
three species. For the three northern species of E. fulvus, 
E. albifrons and E. sanfordi we included an additional 
model of only panmixia of E. albifrons and E. sanfordi and 
only gene flow among the latter two species, as suggested 
by the results of the species tree (Figure 3) and the nuclear 
genetic structure (Figure 4). 
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